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Abstract 

Solving large-scale optimization on-the-fly is often a difficult task 
for real-time computer graphics applications. To tackle this chal¬ 
lenge, model reduction is a well-adopted technique. Despite its use¬ 
fulness, model reduction often requires a handcrafted subspace that 
spans a domain that hypothetically embodies desirable solutions. 
For many applications, obtaining such subspaces case-by-case ei¬ 
ther is impossible or requires extensive human labors, hence does 
not readily have a scalable solution for growing number of tasks. 
We propose linear variational subspace design for large-scale con¬ 
strained quadratic programming, which can be computed automat¬ 
ically without any human interventions. We provide meaningful 
approximation error bound that substantiates the quality of calcu¬ 
lated subspace, and demonstrate its empirical success in interactive 
deformable modeling for triangular and tetrahedral meshes. 

CR Categories: 1.3.5 [Computer Graphics]: Computational Ge¬ 
ometry and Object Modeling 1.3.6 [Computer Graphics]: Method¬ 
ology and Techniques—Interaction techniques 

1 Introduction 

In computer graphics realm, solving optimization with a substan¬ 
tially large amount of variables is often an expensive task. In order 
to speed up the computations, model reduction has been introduced 
as a useful technique, particularly for interactive and real-time ap¬ 
plications. In solving a large-scale optimization problem, it typi¬ 
cally assumes that a desired solution approximately lies in a man¬ 
ifold of much lower dimension that is independent of the variable 
size. Therefore, it is possible to cut down calculations to a compu¬ 
tationally practical level by only exploring variability (i.e., differ¬ 
ent solutions subject to different constraints) in a suitably chosen 
low-order space, meanwhile, attempting to produce visually con¬ 
vincing results just-in-time. In this paper, we re-examine model re¬ 
duction techniques for quadratic optimization with uncertain linear 
constraints, which has been widely used in interactively modeling 
deformable surfaces and solids. 

Modeling deformable meshes has been an established topic in com¬ 
puter graphics for years [Sorkine et al. 2004; Yu et al. 2004]. 
Mesh deformation of high quality is accessible via off-line solv¬ 
ing a large-scale optimization whose variables are in complexity 
of mesh nodes. A studio work-flow in mesh deformable model¬ 
ing often involves trial-and-error loops: an artist tries different sets 
of constraints and explores for desirable poses. In such processes, 
an interactive technique helps to save the computation time where 
approximate solutions are firstly displayed for the purpose of guid¬ 
ance before a final solution is calculated and exported. Neverthe¬ 
less, interactive techniques related to real-time mesh modeling has 
been less successful than their off-line siblings till today. Exist¬ 
ing work based on model reduction often requires a high quality 
subspace as the input, which typically demands human interven¬ 
tions in constructing them. Exemplars include cage-based defor¬ 
mations [Huang et al. 2006; Ben-Chen et al. 2009], LBS with addi¬ 
tional weights [Jacobson et al. 2012], LBS with skeletons [Shi et al. 
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2007], and LBS with point/region handles [Au et al. 2007; Sumner 
et al. 2007]. The time spent on constructing such reduced models 
is as much as, if not more than, that spent on on-site modeling. In 
industrial deployments, companies have to hire many artists with 
expertise skills for rigging a large set of models before those mod¬ 
els are used in productions. This poses the necessity for a fully 
automatic subspace generation method. This problems have re¬ 
ceived attentions in the past. Eor example, data-driven methods 
have been developed for deformable meshes, where a learning al¬ 
gorithm tries to capture the characteristics of deformable mesh se¬ 
quences and applies to a different task [Sumner et al. 2005; Der 
et al. 2006]. However, they still struggle to face two challenges: 
1) Scalability: Like approaches relying on human inputs, obtaining 
a deformable sequence of scanned meshes can also be expensive. 
No words to say if we want to build a deformable mesh database 
containing large number of models with heterogeneous shapes. 2) 
Applicability: Many models of complex geometries or topologies 
are relatively difficult to rig, and there are no easy ways to build a 
set of controllers with skinning weights to produce desirable defor¬ 
mations. Though we see there have been several workarounds for 
a domain specific mesh sets, such as faces and clothes, an automat¬ 
ically computed subspace for arbitrary meshes, which is cheaply 
obtained, still can be beneficial, if not all, for fast prototyping or 
exploratory purposes: the set of constraints chosen on-site is ex¬ 
ported for computing a deformation with full quality in the off-line 
stage. 

In this paper, we introduce an automatic and principled way to cre¬ 
ate reduced models, which might be applied to other computation¬ 
ally intensive optimization scenarios other than mesh deformation. 
Our main idea is very simple: in solving a constrained quadratic 
programming, we observe that Karush-Kuhn-Tucker (KKT) con¬ 
dition implicitly defines an effective sub space that can be directly 
reused for on-site subsequent optimization. We name this linear 
variational subspace (for short, variational subspace). Our contri¬ 
bution is to theoretically study the approximation error bound of 
variational sub space and to empirically validate its success in in¬ 
teractive mesh modeling. The deformation framework is similar 
to one used in [Jacobson et al. 2012].^^ We further examine the 
deformation property of our proposed method, and compare with 
physically based deformation [Ciarlet 2000; Grinspun et al. 2003; 
Botsch et al. 2006; Sorkine and Alexa 2007; Chao et al. 2010] and 
conformal deformation [Paries et al. 2007; Crane et al. 2011]. 

2 Mathematical Background 

Consider minimizing a quadratic function f{X;q) = 
HX — X subject to linear constraints X = b, 
where X G is an overly high dimensional solution, H is 

Tn independent work reported in a recent preprint [Wang et al. 2015], 
Wang et al. also propose a mesh deformation framework based on linear 
variational subspace similar to ours with the difference that we in addition 
use linear variational subspace to model rotation errors in reduced-ARAP 
framework. Our deformation can be similar to theirs [Wang et al. 2015], if 
regularized coefficient a is set to a large value. Therefore, the contribution 
of our paper excluding the empirical efforts is the approximation theory for 
linear variational subspace. 

^Implementations and demos: https : / /github . com/bobye 




Figure 1: Variational sub space provides robust and high quality 
deformation results regarding arbitrary constraints at runtime. This 
figure shows tree and fertility as well as their deformed versions. 



Figure 2: Artistic freedom is important in deformable mesh mod¬ 
eling because many artists are interested in authoring creative 
editing. Rig or cage based deformation lacks the richness of de¬ 
formable variants. 


a semi-definite positive matrix of size n x n, q G ^ is a 
well-conditioned matrix of size n x m, and b G R"^. Typically, 
m <^n. Without loss of generality, we can write 

minx [\/2)X'^HX -q^X 
s.t. A'^X = h. 

Instead of solving the optimization with a single setup, we con¬ 
sider a set of them with a prescribed fixed H, and varying A, b 
and q under certain conditions. The “demand” of this configura¬ 
tion is defined to be a particular choice of A, b and q. Different 
choices usually result in different optimum solutions. When n is 
relatively small, efficiently solving for unreduced solutions belongs 
to the family, so called multi-parametric quadratic programming, 
or mp-QP [Bemporad et al. 2OO2][T0ndel et al. 2003]. We instead 
approach to tackle the same setting with a large n by exploring ap¬ 
proximate solutions in a carefully chosen low-order space. 

We model the “demand’s by assuming each column of Anxm 
is selected from a low-order linear space Cnxd, namely A = 
CnxdAc G Span(C) for some Ac, and q is again selected from an¬ 
other low-order linear space Span(T)) such that that q = DnxkY 
for some Y, where Ac is a matrix of size d x m, V is a. vector of 
size k X 1. Here d and k is the dimension of reduced subspace C 
and D articulating to what A and q belong, respectively. Instead 
of pursuing a direct reduction in domain of solution X, we analyze 
the reducibility of “demand” parameters A and q by constructing 
reduced space C and D. Specifications of the on-site parameters 
Ac, b and Y turn out to be the realization of “demand”s. We can 
rewrite (1) as 

minx.z {1/2)X^ HX-Y^Dl^^X 
s.t. = AtZ = b. 

Optimization (2) can be decomposed into an equivalent two-stage 
formulation, i.e.. 


min { min f(X;DY)}. (3) 

ATz=hCTx = Z 


Karush-Kuhn-Tucker (KKT) condition yields that the optimum 
point X*{Z, Y;C,D) for mm(jTx=z fi^'^DY) should satisfy 
linear equations 


H C \ f X^ \ f DY \ 

o)\a)~\ Zj- 

where A is a Lagrange multiplier. Therefore X"^ {Z^Y]C, D) is 
affine in terms of on-site parameters Z and Y, i.e.. 


( 5 ) 


where N{;C) and U{;C)D can be computed before Ac, b and 
Y are observed in solving the second stage of (3): From Eq. (4), 
each column of iV(; C) is computed through a preconditioned lin¬ 
ear direct solver by setting Z as the corresponding column of Idxd 
with y = 0; And similarly, each column of U{]C)D is linearly 
solved by setting DY to be the corresponding column of Dnxk 
with Z = 0. Remark it is particularly required that C and D has to 
be full-rank and well-conditioned (as will be specified later). 

Once we have a subspace design (5), for arbitrary “demand” A, 
b and q, we can immediately solve for an approximate solution 
-A*(Zmin, F^min) via Substituting (5) into the original formulation 
(1). We clarify that if our assumption is held, i.e., A = CAc and 
q = DY for some Ac and Y, the approximate solution is identical 
to the exact optimal X^in- Furthermore, the approximated solution 
can work with hard constraints in the online solvers. 

In next, we demonstrate its effectiveness by implementing an in¬ 
teractive mesh deformation method based on the model reduction 
framework we proposed. In the end, we will return to the theoreti¬ 
cal aspects, and derive an error bound for the approximate solution 
with respect to the use of C and D. 

3 Interactive Mesh Deformation 

In this section, we start from the point that is familiar to the graph¬ 
ics audiences, and proceed to the practice of our reduced model, 
where we mainly focus on deriving the correct formulation for em¬ 
ploying variational subspace. Experimental results are provided in 
the end. In order for practitioners to reproduce our framework, we 
describe the details of our implementation in Appendix 5. We re¬ 
mark that, subspace techniques described in this section has been 
standardized as described in [Jacobson et al. 2012]. The main dif¬ 
ference is to replace the original linear subspace of a skinning mesh 
with the variational subspace described in our paper. Our varia¬ 
tional subspace techniques extends the fast deformable framework 
as proposed in [Jacobson et al. 2012] to meshes whose skinning 
is not available or impossible, such as those of complex typologies. 
There are, however, good reasons to work with linear blending skin¬ 
ning, for example, it is often possible for artists to directly edit the 
weights painted on a skinned mesh. 

Notations. Denote by vi,..., G R^ the rest-pose vertex posi¬ 
tions of input mesh AT, and denote the deformed vertex positions 
by vi,...,v; € M®. 

Use bold lowercase letters to denote single vertex v G R^ and 3x3 
transformation matrix r, and bold uppercase letters V and R to de¬ 
note arrays of them. We use uppercase normal font letters to denote 
general matrices and vectors (one column matrices) and lowercase 
normal letters to denote scalars. We may or may not specify the 


W* (Z, y ; C, D) = N{; C)Z + U{; C)DY , 







Figure 3: Results subject to different regularization coefficients. 
Deformed solid cylinder upon three point constraints. From left to 
right: a = 0.01, 0.05, 0.1,1,10. 


dimensions of matrices explicitly in the subscripts, hence Mnxn 
and M are the same. For some other cases, subscripts are enumer¬ 
ators or instance indicators. We use superscripts with braces for 
enumerators for matrices, e.g., 

Use II'll to denote Frobenius norm of matrices (vectors), 11-112 to 
denote L 2 norm of matrices, and || ^|| = ^tr(Z^MZ) to denote 
Mahalanobis norm with semi-positive definite matrix M. 



Figure 4: Results with varied number of rotational proxies. From 
Left to right: d = 5, 9,17.‘ for each, deformed result with 32 rota¬ 
tional proxies is shown with dark contours. 


cuts surface/solid mesh into d patches. We revise the original en¬ 
ergy by assuming 

Yk ~ + Qfe , (8) 

where G S'0(3) denotes i/c-th patch-wise frame rotation of the 
cluster that the vertex k belongs, and 



/ ^0 

-qs 

q2 

qfc = 

qs 

qo 

-qi 


V-^2 

qi 

qo 




(9) 


Denote dot product of matrices by o, and Kronecker product of 
matrices by (g). Let I be the identity matrix, 0 be the zero matrix 
and 1 be a matrix of all ones. 

3.1 Variational Reduced Deformable Model 


It should be noted that Laplacian surface editing [Sorkine et al. 
2004] utilizes to approximate the rotation matrix, whereas we 
use q/e to approximate the difference of two rotation matrices 
rfc — . It leads to a different energy by appending L2 penalties 

subject to the regulators q^: 


ARAP energy. In recent development of nonlinear deformation 
energy, the As-Rigid-As-Possible energy [Sorkine and Alexa 2007; 
Xu et al. 2007; Chao et al. 2010] is welcomed in many related 
works, in which they represent deformations by local frame trans¬ 
formation. The objective energy function under this representation 
is quadratic in terms of variables: vertices and transformation matri¬ 
ces with orthogonality constraints. This family of energy functions 
can be written as 

^(v',r) = ^E E 


£'{V\ Q, S) = £iV\ S+Q)+a E 11^^= f+/5 y] f , 

k=l k=l 

( 10 ) 

where Ok denotes the element-wise area/volume, a denotes the 
penalty coefficient of overall spatial distortion, and /3 denotes the 
additional penalty coefficient of surface normal distortion (if appli¬ 
cable). a and /3 are empirically chosen. (See Fig. 3 for deforma¬ 
tions subject to different penalty coefficients a + /3). One potential 
issue is that using this penalty may incur surface folds when shapes 
are bended at large angles. To counteract such effects, we option¬ 
ally use an extra regularization term appended to S' (V', Q, S) pe¬ 
nalizing the moving frame differentials [Lipman et al. 2007], i.e.. 


where Gk are their corresponding sets of edges (see Figure 5 
of [Jacobson et al. 2012]), ajk G M are typically the cotangent 
weights [Chao et al. 2010], and Vk G SO{3) denotes the local 
frame rotations. By separating quadratic terms and linear terms 
w.r.t. Vi, and vectorizing (vi)7=i and {rk)k=i to Vsnxi andi^grxi 
respectively, ARAP energy can be further expressed as 

E(y', R) = - R^KV' + constant, (7) 

where r^r/e = hxs (see [Jacobson et al. 2012] for more details). 


where FL is the set of neighboring local frames and Okj = {ak + 
aj)/2. The two bending cylinder examples in Fig. 8 are produced 
by penalizing the moving frame differentials. 

Let Sgdxi and Q^rxi be the vectorization of (si) and (q^) respec¬ 
tively. S' is quadratic in terms of V' and Q, and its partial gradient 
w.r.t. V' and Q is again linear in terms of S. Hence again we can 
write S' as 


Rotational proxies. By observation, minimizing ARAP energy 
involves solving R, which is in complexity of mesh geometries. 
We modify the original ARAP energy to a piece-wise linear form, 
which relieves the high non-linearity of optimization, but simulta¬ 
neously increases the complexity by introducing linearization vari¬ 
ables. 

We divide r local rotations into d rotational clusters spatially, which 
is an over-segmentation of input meshes. Rotations within each 
patch segment are desired to be similar in deformations. Em¬ 
pirically, we found that a simple k-means clustering on weighted 
Laplacian-Beltrami eigen-vectors fits well with our scheme, which 


E'(y',Q,S) = ^[V-,Qf L[V-,Q]-S'^M[V-,Q]+S'^NS+constmt, 

where S are the rotational proxies of our model, N / 0 iff the extra 
regularization (11) is present. 

There is an interesting discussion about the difference between 
S and S', because S' includes near-isotropic scaling which has 
arguable values over distortion in only one direction for artistic 
modeling purpose in case the desired deformation is far from the 
isometry [Sorkine et al. 2004; Lipman et al. 2008]. (See Fig. 9 for 
comparison with the ARAP energy.) 



Linear proxies. Besides rotational proxies, we add 3m linear prox¬ 
ies via pseudo-spatial linear constraints, 


( 12 ) 

where X are the linear proxies of our model. 

Intuitively, W = Wmxn ^ hxs spans a finite dimensional linear 
space to approach the uncertainty set of onsite constraints provided 
by users. Its choice refiects how we reduce the dimension of antici¬ 
pated constraints, as suggested by the use of variational subspace, A 
simple choice is a sparse sampling of m vertices (shown as Fig. 5), 
i.e (under a permutation) 


sample at singe vertex i 

WsnXSm = {...;0,...,0, l,0,...,0;...}nxm 0-^3 X3, 

and an alternative one is to utilize m vertices groups via clustering, 
i.e., (under a permutation) 


vertices group j of size A/* 


W3nX3m = ' • 7 1 5 0 ? • 


) Fax 


j-th rows 


To this point, technically contrast with our approach, standard 
model reduction technique employ a strategy that vertices are ex¬ 
plicitly represented in low-order by Va'^xi = Ffsnxsm^smxi- In 
order to compute a reasonable subspace, different smoothness crite¬ 
rion are exposed on computing K, such as heat equilibrium [Bar an 
and Popovic 2007], exponential propagating weights [Jacobson 
et al. 2012], biharmonic smoothness[Jacobson et al. 2011]. We 
instead reduce the dimension of constraints, and the subspace are 
then automatically solved accordingly. 


Variational subspace. With context of approximated energy E', 
we are to solve the linear variational problem so as to derive a re¬ 
duced representation of V' in terms of proxies S and X, i.e., 

min E'iV',Q,S) 

(13) 

S.t. IV3mX3nV'=W. 


By KKT condition introducing Lagrange multipliers A, we have 
a set of linear equations in respect of V', Q, A, which can be ex¬ 
pressed as matrix form 

{w = (14) 

This then implicitly establishes a linear map 

[V'-Q] = NwX + UwS, (15) 

where each column of matrices Nw^Uw can be pre-computed 
by a sparse linear solver with a single preconditioning (LU or 
Cholesky), subject to each single variable in vector X and S. 
Solving for Nw and Uw only need one time computation in the 
offline stage. 


Sub-manifold integration. Provided variational subspace, X and 
S span a sub-manifold of deformations. We then restrict our scope 
to determine reduced variables X and S. We employ a routine 
similar to alternating least square [Sorkine and Alexa 2007], where 
we alternatively update X and S via two phases. 


Phase 1: provided solve for X^^\ 

By substituting (15) into approximated ARAP energy (12), we de¬ 
rive a reduced ARAP energy as 

E"{S, X) = ^X'^LX - S'^MX + constant. (16) 

where L = {NwYLNw and M = MNw + {Uw^LNw- With 
onset hard constraints WeqV' = Peq specified by the user (where 
Weq are positional constraints and Peq are their values), we are then 
to solve for linear proxies X 

min E”(S^^\X) 

^ (17) 

s.t. VeqX = Peq - 

where Aeq = Weq Aw, Peq = WeqPw. Hard constraints are the 
default setting of our framework. 

Alternatively, we can pose on-site soft constraints as 

rmn E" (5«,V:)+5||iVeqV: + (7eq5« -Peq|p, (18) 

where (5 > 0 is adjusted interactively by user to match the desired 
effects. We input Weq, Peq and solve the integrated reduced model 
interactively. 

Remark that because optimization problems (equations (17) and 
(18)) are again linear variational, it can be efficiently solved by a 
standard dense linear solver: (1) pre-computing LU factorization 
of matrix (not related to S) at the stage to specify constraint 
handlers Weq, and (2) backward substitution on the fiy at the stage 
to drag/rotate handler. 


Phase 2: provided X^'"^ and S^''\ compute Rather than 

minimizing the reduced energy functional E" (shown in Eq. (16)) 
in terms of S, we instead want rotational clusters to adapt for the 
existing deformation. Letting = NwX^'^^ + UwS^^\ 

we fit a patch-wise local frame of rotational clusters subject to de¬ 
formed mesh by dumping relations and their penalties 
a = = 7 = 0, and optimizing a simplified energy 8 (, S) = 
E'(V'^^^ 0,6'), which is equivalent to 

max 0] = S^iMNX^^ + MuS^^) 

s L . J V 7 

s.t. Si G 60(3), z = l,...,(i, 

where Mn and Mu are pre-computed. It is well known that 
those rotation fittings can be solved in parallel via singular value 
decomposition of each gradient block of Si. For 3x3 matrix, 
we employ the optimized SVD routines by McAdams and col¬ 
leagues [McAdams et al. 2011] that avoid reflection, i.e., guarantee 
the orientation. 


3.2 Algorithm overview 

We review our previous mathematical formulations, and summarize 
our algorithm into three stages (see also Fig. 5): 

Pre-compute. The user loads initial mesh model Al, linear proxies 
W, rotational proxies (ik), and affine controllers (if applicable). 
Our algorithm constructs a sparse linear system to solve for 
variational subspace Nw and Uw, and then pre-computes L, M, 
Mn and Mu- 
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Figure 5: Interactive mesh modeling framework of our approach. From left to right: (1) Input red demon model V; (2) Pre-process to set 
pseudo-constraint points as linear proxies and near-rigid parts as rotational proxies; (3) Compute variational subspace Nw, Uw offline, 
and load into display device. (4) Prepare subspace integration L, M, Mn and Mu- (5) At run-time, given on-site user demand Weq, Peq, 
solve for reduced variables X, S and upload to display device (bandwidth saving); (6) Display deformed mesh and feedback. 



Figure 6: Our reduced model preserves the nature of different func¬ 
tions. Left: volumetric deformed mesh; Right: deformed surface 
(self-occlusion possible) under same constraints. 


Prepare on-site constraints. When above pre-computed matrices 
are present, the user can only freely specify the intended constraint 
handler on-site. They are in the form of Weq. Our algorithm then 
proceeds to compute iVeq and f/eq, and pre-factorize the linear 
system (see equation (17) or (18)). If a user introduces a brand new 
set of constraints on-site, this stage will be re-computed within tens 
of milliseconds. The timing regard to different settings has been 
reported in Table 1 column “OP”. 

Deform on the fly. Our algorithm allows the user to deform 
meshes on the fly, which means the user can view the deformation 
results instantly by controlling constraint handlers. For each frame, 
our model takes in positional constraints Peq, calls an alternating 
routine (with global rotation adaption) interactively to solve for 
proxy variables X and S, and reconstructs and displays the 
deformed mesh. To guarantee real-time performance, we used a 
flxed number of iterations per frame. By initializing an alternating 
routine with the previous frame proxies, we do not observe any 
disturbing artifacts even when using only 8 iterations. 


3.3 Results and Discussion 

We implement our framework for deformable mesh modeling and 
demonstrate our results by examples which include standard de¬ 
formation suites introduced in [Botsch and Sorkine 2008]. Re¬ 
sults of our approach on a set of typical test meshes are shown in 
Fig. 8). The results shown can be compared with results of high 
quality methods without model reduction, including PriMo[Botsch 



Figure 7: Our revised energy does not reveal linearization artifacts 
when the deformation is isometric (middle). It also favors isotropic 
scaling when the model is stretched or shrunk (right). The original 
model is shown in the left. 


et al. 2006], also shown in [Botsch and Sorkine 2008]. Besides, 
we also demonstrate the strength of our method in conformal set¬ 
ting, where we conflgure scaling factors in our modeling frame¬ 
work (see Fig. 10). In our experiments, the modeling framework 
runs robustly on various models, for different types of transfor¬ 
mation, such as small and large rotations, twisting, bending, and 
more as shown in Fig. 2. It works reasonably naturally on both sur¬ 
face and solid meshes, in which user’s choice of energy controls 
the desired behaviors (see Fig. 6). It also accommodates different 
hyper-parameter setting, such as the number and type of proxies, to 
produce predictable and reasonable results (see Fig. 4). 

Based on our CPU implementation. We report timing of our algo¬ 
rithm working on different models presented in our paper in Ta¬ 
ble 1. All timing results are generated on an Intel Core™ i7- 
2670QM 2.2GHz x8 processor with 12 GB RAM. It has been 
shown that the time used in reduced model iteration is not related to 
the geometric complexity, and the overall computation per frame is 
magnitude faster than that as reported in [Hildebrandt et al. 2011]. 
The computational framework used in our paper is almost as same 
as the one used in [Jacobson et al. 2012], thus the performances are 
comparable. It is also shown that the process of mesh reconstruc¬ 
tion, which is a matrix-vector product (see Eq. (15) and column 
“Df.” of Table 1), is the bottle-neck of overall computation, yet it is 
embarrassingly parallel. 

Comparing to other cell-based model reduction methods [Sumner 
et al. 2007; Botsch et al. 2007], our approach utilizes a much 
smaller number of reduced variables. Typically, we adopt no more 
than 35 linear proxies, and no more than 60 rotational proxies. 



















Input Model 

Proxies 

Runtime 

Pre-computation 


Model 

Vert. 

Type 

Linear 

Rot. 

1 Iter. 

(ms) 

DL 

(ms) 

Total 

(ms) 

Subspace 

(s/GB) 

(ms) 

Fig. 

Cylinder 

5k 

Tri. 

33 

12 

50 

1.4 

2 

14 (0.3) 

14 

8 

Cactus 

5k 

Tri. 

33 

27 

90 

2 

3 

18 (0.4) 

16 

8 

Bar 

6k 

Tri. 

33 

52 

93 

3.8 

4.8 

36 (.5) 

20 

8 

Bumpy Plane 

40k 

Tri. 

33 

27 

85 

14 

15 

200 (2.7) 

38 

8 

Plate Box 

4k 

Tet. 

25 

25 

62 

2 

2.7 

30(.4) 

11 

6 

Solid Cylinder 

8k 

Tet. 

33 

52 

no 

4 

5 

90 (.8) 

22 

3 

Tree 

3.6k 

Tri. 

60 

60 

137 

3.5 

5 

34(.4) 

10 

1 

Fertility 

25k 

Tet. 

29 

26 

78 

5.2 

6 

148 (2.4) 

28 

1 

Dinosaur 

21k 

Tri. 

46 

34 

108 

12 

14 

115(1.7) 

24 

7 

Dragon 

53k 

Tri. 

20 

20 

55 

14 

15 

198 (2.7) 

64 

10 

Red Demon 

80k 

Tri. 

30 

28 

90 

35 

36 

498 (5.8) 

106 

5 


Table 1: Model statistics and serial performance on a HP laptop with an Intel i7 2.20GHz x8 Processor From left to right: number of 
vertices, type of mesh, number of linear proxies, number of rotational proxies, time in p seconds for one iteration, time in milliseconds for 
mapping reduced solution info full space (Intel MKL on CPU performance), time in milliseconds for full optimization per frame, time in 
seconds (and memory in GBytes) for pre-computation of subspace, time in milliseconds for computation of on-site preconditioner, figure that 
shows the configuration. 



Figure 8: Approximation quality of our method in all images is demonstrated on the test suite of models introduced in [Botsch and Sorkine 
2008]; From left to right: Bumpy Plane, Cylinder, Cactus, Bar. Multiple types of deformations, including bending, shifting and twisting, are 
tested. 


Besides, the configuration of proxies gives user the freedom to 
design his/her own needs in modeling a particular mesh. Instead 
restricting variability in modes and modal derivatives space [Hilde- 
brandt et al. 2011], artist, based on his intentions, can cut shape into 
near rigid parts (each for a rotational proxy), and specify pseudo 
constrain locations as linear proxies. Fig. 5 demonstrates a model¬ 
ing scenario where artist intended to adjust mouth, nose and eyes 
on a face model: semantic parts are in first annotated, variational 
reduced deformable model is then pre-computed for on-line editing. 


4 Theory of Variational Subspace 

The notations used follows the preliminary setup in section 2. 

4.1 Concept of Approximation 

Definition 4.1 (Variational Subspace). Given a quadratic program¬ 
ming problem in the form of Eq. (1), choose C and D for the prob¬ 
lem in the two-stage form as Eq. (3). The solution X*{Z, Y;C,D) 
for the first stage problem is hence given by Eq. (5). The subspace 
spanned by columns of N and U-D are called variational subspace. 
Proposition 4.1. N is a matrix of size n x d and U ■ D is a 
matrix of size n x k. Their columns span a linear subspace 
Span(f/)+Span(f/ -D) where the reduced solution belongs. Those 
columns are computed by solving a variational formulation pro¬ 
vided by Eq. (4). 

In Eq. 1, iT is only guaranteed to be positive semi-definite. We have 


Cholesky decomposition H = L, where L is upper triangu¬ 
lar with non-negative diagonal entries [Golub and Van Loan 2012]. 
Denote the pseudo-inverse of L by L^, then for any X G 'MF, we 
have a two-part orthogonal decomposition X = X X, where 

X = L+LX, X = {I- L+L)X. (20) 

Proposition 4.2. With the two-part decomposition Eq. (20), we 
have 

LX = LX, LX = 0, X'^X = 0. (21) 

In addition, if H = LZ L is positive definite, X = 0. 

Moreover, we can rewrite the optimization problem Eq. (1) as 

minx (1/2)(X + X)^L^L(X + X)-g^(X + X) 

s.t. ^^(X + X) = b (22) 

X = L+LX, X = {I- L+L)X. 

Since LX = 0, we simplify above formulation to 

minx f{X) = {l/2)X'^L'^LX-q^X-q^X 

s.t. A^X = b- A'^X (23) 

X = L+LX, X = (7 - L+L)X. 

It is observed that only X appears in second-order term in the 
objective function of Eq (23). Suppose the optimal solution to 
Eq. (1) is Xmin with a two-part decomposition (given by Eq. (20)) 


























Xmin = -^min + -^min, wc then considcr the following companion 
optimization problem 

min^ /(X) = (1/2)X^L'^LX - X - q^Xr^in 

s.t. A^X = b- A^Xmin (24) 

(/ - L+L)X = 0. 

Remark —g^Xmin which appears in objective function of Eq. (24) 
is a constant. We see Eq. (24) can be equivalently solved in two 
steps: In the first step, we solve the following problem 

min^ /(X) + q^^X^in = (1/2)X'^L'^LX - g^X 

s.t. X = b - Xmin- ' 

And in the second stage, we project the solution to L'^LXmin, 
where Xmin is the solution to Eq. (25). 

Theorem 4.3. Suppose Xmin is the unique solution to Eq. (23), 
Xmin is the unique solution to Eq. (24), and Xmin is the unique so¬ 
lution to Eq. (25), then we have X^^^ = L^LXmin = L^LXmin- 

Proof. We observe that L+LXmin satisfies the constraint of 
Eq. (24), and objective functions of Eq. (23) and Eq. (24) coin¬ 
cide, i.e., /(Xmin) = /(L^LXmin). Therefore, as X^i^ is the 
minimizer to Eq. (24), we have 

/(Xmin) > /(X^in)- 

On the other hand, /(X^in) = /(X^in + Xmin). Given Xmin is 
the minimizer to Eq. (23), we also have 

/(X^in) > /(Xmin). 

Hence, the equality holds for /(X^in) = /(Xmin). Given that the 
optimum exists and is unique, we have 

^^min 4” X^min — X^min X^min — T/X^min. 

The proof of X^^^ = L+LXmin is similar: observe that 

/(X^in) = /(L^LXmin) and L^LXmin satisfies the constraint 
ofEq. (24). □ 

Definition 4.2. Two solutions subject to the form Eq. (1) (param¬ 
eterized by A, q, and b) are called quotient equivalent, if they 
share the same companion problem defined by Eq. (25), i.e. their 
h{A^ b,q) — b — A^X^in are the same. This forms group equiva¬ 
lence in the space of solutions. 

Eor example, let H be the Laplacian operator A, the solution would 
minimize the the L2 norm of first-order gradient. In such case, two 
problems are of quotient equivalence if their optimal solutions pre¬ 
serve to an additive constant. We use the distance between two 
solution groups under the quotient equivalence to measure the ap¬ 
proximation error. It is the Mahalanobis distance provided by H, 
i.e. 

dH{x,y) = {x- vYH{x -y) = {Lx - LyY'{Lx - Ly). 

Let LX = X, q = L'^q, A = L'^A, and b = b — A^Xmin, we 
rewrite Eq. (25) (but not equivalent) as 

miny {1/2)X'^X -q'^X 

s.t. A'^X = b (26) 

{I - LL+)X = 0. 


Similar to the treatment of Eq. (24), we can in first solve 

min^ (1/2)X^X - fX 

s.t. A^X = h. ^ ’ 

and then project the solution to LL^Xmin, where Xmin is the so¬ 
lution to Eq. (27). 

Since we are always interested in distance measure dn for different 
solutions, the projection step in solving Eq. (24) is not necessary to 
compute dn- The distance between two solution groups Xi and 
X 2 of the quotient equivalence is therefore the Euclidean distance 
between Xi = LXi and X 2 = LX 2 . 

Similar to the treatment of Eq. (1) and Eq. (27), we can derive the 
two-stage problem from Eq. (2) as 

minj> (l/2)X^X-y^£)^X 

s.t. C'^X = Z (28) 

AlZ = h, 

where D = {L+fD and C = {L+^C. The KKT condition of 
its first-stage problem is given similarly as 



where A is a Lagrange multiplier. We have the following justifica¬ 
tions to only study Eq. (27) and Eq. (28). 

Proposition 4.4. If Xmin is the optimal solution to Eq. (1) and 
^Amin the optimal solution to Eq. (27) with Xmin = (/ - 

L 7/)-A^min and b — b A. .^Amin? y^e have L-A^min — LL -A^min* 
The similar argument also holds for Eq. (2) and Eq. (28). 

Proof. Because L^LXmin is the optimal solution to Eq. (23), we 
have L^LXmin = L^LXmin, where Xmin, as mentioned, is the 
optimal solution to Eq. (25). On the other hand, we know LXmin = 
LL+Xmin. Therefore, we have L(L+LXmin — L+Xmin) = 0 ^ 
LXmin = LZ/~*"Xmin- D 

Proposition 4.5. Suppose X* = NZ -\- UDY is the solution of 
Eq. (4) and X* = NZ -\-UDY is the solution ofEq. (29), we have 
LX* =X* if Z = Z - - L+L)X* and Y ^Y. 

Proof We have (L+)^(L^LX* + CA - DY) = 0 ^ LX* + 
CA - L)y = 0 and C'^L^LX* = Z - {I - L+L)X* ^ 
LX* = Z. Therefore, LX* satisfies Eq. (29). □ 

Definition 4.3 (Variational Subspace Under Quotient Equivalence). 
Given X* = NZ UDY to be the solution to Eq. (29), the 
columns of N and tj • D span the variational subspace Span(iV) + 
Span(f/ • D) for solutions of optimization problem in the form of 
Eq. (27). 

The main problem is thus revealed: how close are the exact so¬ 
lution of Eq. (27) and subspace solution restricted in a variational 
subspace defined by Def. 4.3, where the closeness is measured by 
the Mahalanobis distance dn provided by if = L^ L. 



4.2 Bound of Approximation Error 

In this section, we provide the proof that the approximation error of 
model reduction by variational subspace can be bounded in dn • 
Proposition 4.6 (Exact Solution). Assuming Eq. (27) has a unique 
solution which has finite optimum, such that A is a full-rank matrix. 


the solution is 

Xmin = (/ - AA+)q + {A+fh, (30) 

where — {A^A)~^ is the pseudo-inverse of A. 

Proof The KKT condition of Eq. (27) indicates Xmin — q — 
AAmin and A^Xmin = b. This leads to A^{q — AAmin) = b. 
Since A is full-rank, A A is invertible. Hence we have 

Arain = {A^A)-\A^q-b). (31) 

Plug-in Xmin = q - AAmin yields Eq. (30). □ 

Similarly, we also have 

Proposition 4.7. Suppose X*(Z, Y) = NZ UDY is the solu¬ 
tion of Eq. (29), then 

N = (6+)^, UD = {I- CC+)D, (32) 

where C+ = is pseudo-inverse of C, and 


U = I — CC^ is the orthogonal projector onto the kernel of 
[Golub and Van Loan 2012]. 

Here we remark that in order for any subspace solution X (Z,Y) 
to have a unique low-dimensional coordinate {Z,Y). We should re¬ 
quire C and D to be linearly independent. This equivalently means 
UD is full rank. 

Theorem 4.8 (Projection on Variational Subspace). Assume 
columns ofC and D are linearly independent. Given any X G 
its closest point (under Euclidean distance) in a variational sub¬ 
space X*(Z, Y) given by Eq. (32) is 

Y = (UD)+X, Z = C'^X, (33) 

and the closest point is 

X* = {UD{UDy + CC+)X, (34) 

where {&£))+ = (D'^UD)-^ D'^U’^ andtj = I - CC+. 

Proof. If UDv = 0 for some v e we have Dv + Cu = 0, 
where u = C'^v. Since C and D are linearly independent, we 
have u = 0 and u = 0. Therefore, UD is full-rank. Eurthermore, 

D^UD = D U UD is invertible. The closest point to X is to 
minimize 

min||X*(Z,y)-X||^ 

Z,Y 

whose partial gradient against Z and Y should be zero, i.e., 
N'^(NZ + UDY - X) = 0 

ifu'^{NZ + UDY - X) = 0. 


Notice that N'^tj = 0, iV^X = (C'^C')“\ (7^(7 = U. Above 
equalities can be simplified to 

(c^cy^z = N^x, b^uDY = b^u^x. 

Above equalities can be solved as 

Z = {C'^C)N'^X = c'^x, Y = {D'^UDy^D'^U'^X. □ 

Next, we are to derive the analytic subspace solution. 

Proposition 4.9 (Variational Subspace Solution). Let I = 
Ub(Ub)~^ + CC~^ be orthogonal projector onto the subspace 
Span(C) + Spsin(b)([Yanai et al. 2011], page 45), where U = 
/ — CC^ the orthogonal projector onto the kernel ofC^. Assum¬ 
ing A is full-rank, columns of C and D are linearly independent, 
and (A^IA)~^exists, the variational subspace solution to Eq. (27) 
is 

= (/ - Aj,A+) q + {A+yi (35) 

where Ap = IA and Ap = ^Ap Ap^ Ap. Note 7 — ApAp is 

the projection matrix restricted in subspace Span((7) + Span(l)) 
that map onto the kernel of Ay 

Proof. Eirst, we have / is symmetric, and II — I. Plug variational 
subspace X*(Z,Y) into Eq. (27). Erom the KKT condition, we 
have 

( + UDYrnin + iAmin “ «] = 0 

and 

A'^ [NZ^,n + UDYrpf = b. (36) 

Similar to the derivation in Theorem 4.8, the former equality of 
KKT condition yields 

= {UDy{q - iAX„), = C'^iq - AA^J, 

and 

XYn = i{q-AAXy. (37) 

^ ^ 

Let Xp = q — AAy^, and combine Eq. (37) with Eq. (36), we 
have A^i(q — Akyy — b. It gives 

A:,in = {A^iA)-\A^iq - b). (38) 

Plug Eq. (38) back to Xy^ (Eq. (37)) yields the subspace solution 
Eq. (35). □ 

We are now ready to introduce the main result. Let || • || be the 

induced L 2 matrix norm, which is its largest singular value. 
Theorem 4.10 (Approximation Error Bound of Variational Sub¬ 
space Solution). Given the demand matrix C and b forming the 
subspace Span(C) + Span(i9), where U — I — CC^, and 
I — Ub(Ub)^ CC^. The error between reduced solution 
-Amin io Eq. (35) and exact solution Xmin to Eq. (30) has a fol¬ 
lowing upper bound: Assuming ||/ — < p < 1 and 

cond(A) = ||A||||A+|| < CU < +(X) for any A in the scope of 
optimization Eq. (27), there exists constants > 0 and ^2 > 0, 
such that 

Ifc-x^inll < \\iq-q\\ + x(bA,A+-y,y ||7i-i||, 

(39) 


and 



for any q, b, and C, D, A, where 

A {b,q,A+-fiufi2)=Mb\\ ■ ||i+||^ + /32II4II • ll^+ll > 0. 

In particular, if q = DY and A = CAc for some Y and Ac, then 
it must have Iq = q and IA = A, thus we know = ^min- 

Proof See Appendix 5 □ 

Theorem 4.10 bounds the approximation error between reduced so¬ 
lution and exact solution by two terms. They are the norm of pro¬ 
jections of q and A onto the intersection of kernel space of and 
(7^. Finally, given the Prop. 4.4, we have 

= ||LL+X:,i„-LL+X^in|| 

^ - - (40) 

< ||LL+||||X*i„-X^in|| 

— l|V^min Xminlli 

where X^in is the solution to Eq. (1) and X^^^ is the corresponding 
variational reduced solution. 

5 Conclusions 

In this paper, we presented variational subspace for reducing cal¬ 
culations in minimizing quadratic functions subject to large-scale 
variables, and integrated it into an interactive modeling framework 
for mesh deformations. Variational subspace is an economical sub¬ 
space driven by reduced constraint demands and optimization con¬ 
texts. Based on it, we implemented an easy-to-use mesh manipu¬ 
lator, which is efficent, robust in quality, intuitive to control, and 
extensible. 
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Proof of Theorem 4.10 


Similarly, we have 

\\A* -AH = mX^A)-'- - (A^A,)-^\A'^i\\ + \\A* - i+/|| 

“ 1-||/-A+/A|| . " 

□ 

Corollary. Since ||7 - i+7i|| = \\A+A - A+iA\\ < 

||4+||||74 — A\\, we have 

\\A UA'^A)-^ - (AlA.yA i^ll < cond(7i)p+|| ||;^_^|| 

II r ^ V p p; J II - 1_||J_^+J^||" 


lli+ 


l-||7-i+7i|| 


Here we prove the main result 

Proof Based on Prop 4.6 and Prof 4.9. We can decompose the 
error term into two parts: 

||x;Lin-Xmin|| < \\{Ay^b-{A+)b\\+\\{i-A^A;)q-ii-AA+)q\\. 

(44) 

On one hand, based on Corollary 5, there exists a constant fSi = 

> 0, such that 

i-p 

\\{A+fb-{At)b\\<\\A+-At\\\\b\\ 

<y\\b\\\\A+f\\iA-A\\ 


We follow the notations used in section 4.2 to prove Theorem 4.10. 
We firstly have the following lemma. 

Lemma .1. Let I = UD(UD)^ + CC~^ and Ap = I A, where 
U = I — CC^. Assume ||/ — < 1 (while we know 

II / — A^/A|| < 1, because I — A^ I A is the projection matrix 
onto kernel({A [U £), C]}^) and the equality holds if and only if 
TSink{A~^[UD,C]) < m), we have 


|i UA^A)-^ - (A^Ayf A^W < cond(i)||7-i+7i|| 
I ' V p p; J II - 1 _ ||7_^+7^|| 


(41) 


where cond(4) = ||4|| || is the condition number of A, and 

'|i+|| ||7- i+7i|| 


U+-Aj\\ < 


1-\\I-A+IA\\ 


+ ||7l+r||77l-7l||. (42) 


On the other hand, from Corollary 5, we also have 

\\ii-A^A;)q-ii-AA+)q\\ 

< II(7 - Ai+)4 - (7 - ii+)7g|| + ||7 - ii+ II ||74 - g|| 

< 11(7 - ApA;)q - i{I - AA+)iq\\ + ||i+7g||||7i - ill 
+ ||7-ii+||||7g-g|| 

^ ii^iiM^iiii^^-^11 + m- ^11, 

(46) 

where ^2 = 1 + — > 0- Combining Eq. (45) and Eq. (46), 

1- p 

Eq. (39) is held. □ 

Implementation Details 


Proof We have the following expansion 

i = i(7 - (i+7i)-ii+ = 

i((7 - i+7i) + (7 - A+iAf + (7 - A+lAf + .. .)A+. 

Therefore, we have 

||i - (i^ip)-'] i^ll 

< ||i||(||7 - i+7i|| +1|7 - A+iAf + 1|7 - A+iAf +.. .)||i+|| 

^ cond(A)||/— A^/A|| 


In companion to our proposed algorithm, other algorithm details 
less relevant to variational subspace is provided in this section, 
which we follow the notations used in section 3. 

Global rotation adaption 

In our implementation, we also introduce global rotation ro to di¬ 
minish the approximation error of loeal rotation matrix ineurred 
by piece-wise linear form (see equation (8)). It is fitted again by 
a single SVD in each frame. Then we update the reduced model 
(Phase 1 and 2) under updated frame coordinates, i.e., multiply the 
inverse rotation to Peq and S. Meanwhile during mesh reeon- 
struction (based on Eq. (15)), we should also multiply rotation ro 
to vertiees V', so as to display them in the original frame. 


( 43 ) 




Figure 9: Left: Original model; Middle: Out-of-shape distortion 
in ARAP surface modeling; Right: Our method deforms rest-pose 
part via shrinking. 



spondence, of which the first phase is identical to former. For the 
second phase, we reformulate as follows. 

Similar to the previous Phase 2, we are again to fit the consistent 
local frame by optimizing the simplified energy 


5(V'(‘\S) = constant - (To (^(g>l 9 xi)) 

+1 E E , 

^=1 (iJ)eQk 

(48) 

where Tgdxi is the vectorization of (t^), ^ = [fi,, f)d] and 
S = T o (g) Igxi). To solve for S, we use two steps: first, 
we fix ^ and optimize T, which is exactly the same as discussed. 
It is required to compute and perform sin¬ 

gular value decomposition. Second, we fix T and compute partial 
gradient w.r.t. ^ 


^ = -(/dxd®lix9)ToM[l/'«;0]+Co^ 

= -(/dxd®lix9)To(MivX« +Mu5'«) + Co^ , 

(49) 

where Cdxi is pre-computed and + MuS^^^) is com- 

puted in the former step. Hence by setting = 0, ^ is com¬ 
puted. Fig. 10 illustrates the difference between deformations with 
conformal factors and without. 


Figure 10: Difference between deformations with conformal fac¬ 
tors (right) and without (middle). Seven constraint points are 
marked as red balls. 


Affine Patches 

In deformable modeling, the user usually would like to constrain 
certain patches on the mesh to be rigid or fixed, or more generally 
affine. Our framework can be accompanied by those requirements 
in pre-computation. Vertex positions V' on the deformed mesh can 
be linearly expressed in terms of deformable vertices Vq and patch- 
wise transformation ti, d^, i.e., (under a permutation) 

V' = [V[);Viti+di;...;V.t. + d,] , (47) 

where Vq are deformable vertices, Vi ,..., Vs are s affine patches 
on the original mesh with prescribed transformation matrices t = 
[ti,..., ts], and displacements d = [di,..., ds]. Under this rep¬ 
resentation, the first stage problem is reformulated accordingly such 
that the variational subspace is solved for variables of the de facto 
control layer [Vo,t,d, Q], instead of for [V',Q] (see equations 
(14) and (15)). For simplicity, each affine patch accompanies a sin¬ 
gle rotational proxy and a single linear proxy. 

To improve the numerical stability in case one would like to con¬ 
strain more than one vertex on a single affine patch (e.g., constrain 
four in rigid motion), we in addition append corresponding linear 
proxies for each variable of t. Therefore, the total degree of linear 
proxies is 3m + 9s. 

Conformal-like Deformations 

We extend our model to conformal-like deformations in this sec¬ 
tion, by introducing scaling factors for each rotational proxy. In¬ 
stead of restricting G 5'0(3), we permit ^ S'0(3), where 

||si ||2 G [1/'?/^,'0], for some constant'0 > 1. 

Thus we can write = 'fiti, where 'ipi G and ti G 

SO (A). The updating routine also contains two phases in corre- 



